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ABSTRACT 


A  general  purpose  computer  program  is  developed  to 
perform  nonlinear  constrained  optimization  of  engineering 
design  problems.  The  program  is  developed  especially  for 
use  on*  microcom pute rs  and  is  called  .1  icrocompute r  Software 
for  Constrained  Optimization  Problems  (*SCO?) .  It  will 
accept  a  nonlinear  objective  function  and  up  to  50 
ineguality  constraint  functions  and  up  to  20  bounded  design 
variables. 

M  SCOP  employs  the  method  of  feasible  directions. 
Although  developed  for  microcomputers,  for  speed  of  develop¬ 
ment,  the  KSCOP  was  implemented  on  an  I3M  3033  using  stan¬ 
dard  basic  language,  Waterloo  BASIC  Version  2.0.  It  is 
directly  transportable  to  a  variety  of  microcomputers. 

Typical  applications  of  3SCOP  program  are  in  the  design 
of  machine  components  and  simple  beam  and  truss  structures. 
Solutions  to  three  sample  problems  are  given. 
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I.  INTRODUCTION 

A.  PURPOSE 

This  thesis  describes  the  development  of  a  microcomputer 
oriented  program  called  MSCOP  (tti crocomp u ter  Software  for 
Constrained  Optimization  Problems)  for  constrained  optimiza¬ 
tion  of  engineering  design  problems.  Problems  which  can  be 
solved  by  the  MSCCF  are  nonlinear  programming  problems 
arising  in  several  areas  of  machine  and  structural  design, 
such  as  the  minimum  weight  design  of  structures  subject  to 
stress  and  displacement  constraints  [Ref.  1]. 

In  recent  years,  several  powerful  general  purpose  opti¬ 
mization  programs  have  become  available  for  engineering 
design  problems,  e.g.,  COPES/CONMIN  [Ref.  2],  and  ADS-1 
[Ref.  3].  These  programs  can  handle  a  wide  range  of  design 
problems  and  contain  a  variety  of  solution  techniques. 
Also,  several  programs  are  available  that  include  optimiza¬ 
tion  in  an  integrated  analysis  /  design  code,  e.g.,  ACCESS, 
ASOP,  EAI,  PAPS,  SAVES,  SPAR,  STARS  and  TSO  [Ref.  4].  Ail 
of  the  above  optimization  programs  are  written  in  7CRTF.AN, 
and  are  built  for  use  on  a  mainframe  computer.  Their  use  can 
be  cumbersome,  especially  for  the  occasional  user.  Since 
many  engineers  are  now  using  microcomputers,  there  is  a  need 
to  develop  an  optimization  program  contained  in  a  microcom¬ 
puter  software  package  for  use  on  microcomputers.  This 
thesis  fills  that  need  by  developing  a  compact  program 
written  in  a  standard  EASIC  language  suitable  for  a  wide 
range  of  microcomputers. 
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B.  I EPLEHSNTATTON 

The  nature  of  an  optimization  program  depends  on  the 
computer  and  programming  method  available.  "he  MSCC?  soft¬ 
ware  is  designed  for  use  on  a  microcomputer .  However,  for 
the  speed  of  development  and  testing,  M5C0P  was  developed  on 
the  I  EM  30  33  computer  at  the  F.  F.  Church  Computer  Center  in 
Naval  Postgraduate  School,  and  was  written  in  WEA5IC 
(Waterloo  Basic)  Version  2.0. 

To  make  sure  that  the  program  is  easily  portable  to  a 
microcomputer,  only  standard  BASIC  commands  and  functions 

are  use).  Tor  example,  TOR  I  =  1  TO  MD3 _  NEX"  I,  G0SU3 

etc.,  were  used.  The  commands  and  functions  not  available 
in  all  variations  of  BASIC  are  avoided,  for  example,  TEN  (A), 
MAT  (A) ,  etc. 

MSCOE  provides  design  engineers  with  a  convenient  tool 
for  optimization  of  engineering  design  problems  with  up  to 
20  bounded  design  variables  and  as  many  as  50  inequality 
constraints. 


C.  GENEEAI  OPTIMIZATION  MODEL 

The  general  optimization  problem  to  be  solved  is  of  the 
form  :  Find  the  set  of  design  variables  _X  that  will 

Minimize  F  (X)  ( 1.  1) 

Subject  to  G  (X)  <0  j  =  1,...,  m  (1.2) 

j 

1  u 

X.<  X.<  X.  i  =  1 - -  r.  (1.3) 

i  -  l  -  i 


where  X  is  referred  to  as  the  vector  of  design  variables. 


variables.  Although  these  bounds  or  "side  constraints" 
could  be  included  in  the  inequality  constraint  set  gi /en  by 
Eq  ( 1 . 2)  ,  it  is  convenient  to  treat  them  separately  because 
of  their  special  structure.  The  objective  function  and 
constraint  functions  nay  be  nonlinear,  explicit  or  implicit 
in  X.  However,  they  must  be  continuous  and  should  have 
continuous  first  derivatives. 

In  general  engineering  optimization  problems,  the  objec¬ 


tive  to  be  minimized  is  usually  the  weight  or  volume  of  a 


structure  being  designed  while  the  constraints  gives  limits 
cn  compressive  stress,  tensile  stress,  Euler  buckling, 
displacement,  frequencies  (eigenvalues),  etc.  [Ref.  5  : 
p.254].  Equality  constraints  are  not  included  because  their 
inclusion  complicates  the  solution  techniques  and  because  in 
engineering  situations,  equality  constraints  are  rare. 

Most  optimization  algorithms  require  that  an  initial 
value  of  design  variables  X°  be  specified.  Beginning  from 
these  starting  values,  the  design  is  iteratively  improved. 
The  iterative  procedure  is  given  by 


q+1  q 
X  =  X  + 


a*  S' 


(1.4) 


where  q  is  the  iteration  number,  S  is  a  search  direction 
vector  in  the  design  space,  and  a*  is  a  scalar  parameter 
which  defines  the  amount  of  change  in  X.  At  iteration  q,  it 
is  desirable  to  determine  a  direction  S  which  will  reduce 
the  objective  function  (usable  direction)  without  violating 
the  constraints  (feasible  direction).  After  determining  the 
search  direction,  the  design  variables,  X,  are  updated  by  Eg 
(1.4)  sc  that  the  minimum  objective  value  is  found  in  this 
direction.  [Ref.  6]. 

Thus,  it  is  seen  that  nonlinear  optimization  algorithms 
for  the  general  optimization  problem  based  on  Eg  (1.4)  can  be 
separated  into  two  parts,  determination  of  search  direction 
and  determination  of  scalar  parameter  a*. 


D.  ORGANIZATION  OF  THIS  THESIS 

This  chapter  has  stated  the  purpose  of  the  thesis  and 
has  put  the  general  concept  of  engineering  optimization  into 
a  preliminary  perspective.  Chapter  2  will  describe  the 
essential  aspects  of  the  optimization  algorithm  used  ir. 
MSCOP  such  as  finding  a  search  direction,  the  one- 
dimensional  search  and  convergence  criteria.  Chapter  3 
describes  program  usage.  In  chapter  4,  there  are  three 
examples  which  are  solved  by  the  MSCOP.  Summary  and  conclu¬ 
sions  are  given  in  chapter  5.  The  program  is  listed  ir.  the 
appen  dix . 


II.  OPTIMIZATION  ALGORITHM 

A.  I NTRCDDCTTON 

There  are  many  optimization  algorithms  for  constrained 
nonlinear  problems  such  as  generalized  reduced  gradient 
method,  feasible  direction  method,  penalty  function  methods. 
Augmented  Lagrangian  multiplier  method,  and  sequential 
linear  programming.  The  feasible  direction  method  is  chosen 
for  development  in  this  thesis  for  three  main  reasons. 
First  it  progresses  rapidly  to  a  near  optimum  design. 
Second  it  only  requires  gradients  of  objective  and 
constraint  functions  that  are  active  at  any  given  ucirt  in 
the  optimization  process  [Hef.  7],  Third,  because  it  main¬ 
tains  a  feasible  design,  engineer  cannot  fail  to  meet  safety 
requirements  as  defined  by  the  contraints.  However,  the 
method  does  have  several  disadvantages  in  that  it  is  prone 
to  "zig-zag"  between  constraint  boundaries  and  that  it  is 
usually  does  not  achieve  a  precise  optimum.  This  method 
solves  the  nonlinear  programming  problem  by  moving  from  a 
feasible  point  (can  be  initially  infeasible)  to  another 
feasible  point  with  an  improved  value  of  the  objective 
value . 

The  following  strategy  is  typical  of  feasible  direction 
method  :  Assuming  that  an  initial  feasible  point  X°  is 
known,  first  find  a  usable-feasible  direction  S.  The  algo¬ 
rithm  for  this  is  similar  to  linear  programming  and  comple¬ 
mentary  pivoting  algorithms.  Having  found  the  search 
direction,  a  move  is  made  in  this  direction  to  update  the  X 
vector  according  to  Eg (1.4)  .  The  scalar  a*  is  found  by  a 
one-dimensicnal  search  to  reduce  the  objective  function  as 
much  as  possible  subject  to  constraints.  That  is  dTN 


feasible  region.  After  updating  the  X°  vector,  the  conver¬ 
gence  test  must  be  performed  in  the  iterative  algorithm.  A 
convergence  criteria  used  in  tnis  is  implementation  are 
described  ir  section  E.  The  general  algorithm  used  in  "SCO? 
is  given  in  Figure  2.1 

B.  SEARCH  DIRECTION 

In  the  feasible  direction  algorithm,  a  usable  -  feasible 
search  direction  S  is  found  which  will  reduce  the  objec¬ 
tive  function  without  violating  any  constraints  for  some 
finite  move.  It  is  assumed  that  at  any  point  in  the  design 
space  (at  any  X)  the  value  of  the  objective  and  constraint 
functions  as  well  as  the  gradients  of  these  functions  with 
respect  to  the  design  variables  can  be  calculated.  Since 
these  gradients  cannot  usually  be  calculated  analytically, 
the  finite  difference  method  Eg(2.1)  is  used  in  HSCOF. 


ar  (X)  i  (x+  ae  )  -  p  (X) 

i 


i 


where  e  is  the  ith  unit  vector 
i 

£  is  a  small  scalar. 

In  MSCOP ,  £  is  0.1"  of  the  ith  design  variable 

In  the  feasible  direction  algorithm,  there  are  usually 
or.e  or  mere  "active"  constraints.  A  constraint  S  (X)  <0  is 

"active"  at  X  if  g (X)  ~  0.  As  shown  in  Figure  2.1,  if  no 
constraints  are  active  the  standard  steepest  descent  direc¬ 
tion  S  =  -  7F  is  used. 


Figure  2.2  Usable- Feasible  Direction. 

Assume  there  are  NAC  active  constraints  at  X..  The  direction 
S  is  "usable"  if  it  reduces  the  objective  function,  i.e., 

VF-S  <  0  (2.2) 

Similarly  the  direction  is  feasible  if  for  a  small  movement 
in  this  direction,  no  constraint  will  be  violated,  i.e., 

VG*S  <  0  (2.3) 

This  is  shown  geometrically  in  Figure  2.2 


2 .  Active  Constraints 

It  is  necessary  to  determine  if  a  constraint  is 
active  or  violated  in  the  feasible  direction  algorithm.  A 
constraint  G  (J)  <0  is  "active"  at  X°  if  G  (X°)  ~  0.  In  order 

to  avoid  the  zigzagging  effect  between  one  or  more 
constraint  boundaries,  a  tolerance  band  about  zero  is  used 
for  determining  whether  or  not_  a  constraint  is  active.  From 

the  engineering  point  of  view,  a  constraint  G  (5)  <  0  is 

active  near  the  boundary  G{X)  =  0  whenever  ACC  <  G  (X)  <  7CC. 

ACC  is  the  active  constraint  criterion  and  VCC  is  the 

violated  constraint  criterion  in  MSCOP.  Assuming  the 

feasible  constraints  are  normalized  so  that  G  (X)  ranges 
between  -1  and  0  for  reasonable  values  of  X,  the  constraint 
G  ( X)  <0  is  considered  active  if  G  (X)  >  -0.1.  The 
constraint  is  considered  to  be  violated  if  G (X )  >  0.004. 

This  is  an  algorithmic  trick  which  improves  efficiency  an  1 
reliability  of  the  algorithm.  However,  since  in  the  one  - 
dimensional  search,  all  interpolations  for  constraint  G (X) 
are  done  for  zeros  of  a  linear  or  quadratic  approx imaticn  to 
G  (X)  in  crder  to  find  a*,  at  the  optimum  the  value  of  active 
constraints  are  very  near  zero,  but  may  be  as  large  as  0.004 
[Ref.  6].  From  an  engineering  point  of  view,  a  0.4  r 
constraint  violation  is  considered  to  be  acceptable. 

3.  Suboptimizat  icn  Problem  and  Push-Off  Factors 

Zoutendijk  [Ref.  8]  has  shown  that  a  usable 
feasible  direction  S  may  be  found  as  follows  : 

Maximize  ^  (2.  4) 

Subject  to  ; 


2F  (X) • S  +  ^  <  0 


a 


5) 


18 


je  J 


(2.6) 


<7  G  (X)  •  S 

S  bounded 


6  £  <  0 
j 


(2.7) 


Where  scalar  ^  is  a  measure  of  the  satisfaction  of  the 
usability  and  feasibility  requirements.  The  scalar  0j  in 
Eg  (2.  6)  is  referred  to  as  the  "push-off'’  factor  which  effec¬ 
tively  pushes  the  search  direction  away  from  the  active 


Figure  2.3  Push-Off  Factor  and  Bounding  of  the  S-Vector. 

constraints.  In  Eg  (2. 6) ,  if  the  push-off  factor  is  zero, 
the  search  direction  is  tangent  to  the  active  constraints, 
and  if  it  is  infinite,  then  the  search  direction  is  tangent 
to  the  objective  function.  It  has  been  found  that  a 


push-off  factor  is  defined  as  follows  gives  good  results 

[ Fef .  5 :  p . 167  ]  : 


0 

3 


g  .  m 

3 


ACC 


where  0.  =  1  . 


(2.3) 


To  avoid  an  unbounded  solution  when  seeking  a  usable 
-  feasible  direction  it  is  necessary  to  impose  bounds  on  the 
search  direction  S.  Cne  method  of  imposing  bounds  on  search 
direction  is  to  impose  bounds  on  the  components  of  S-vector 
of  form  : 


-  1  <  S  •  <  1 


(2.9) 


This  choice  of  bounding  the  S-vector  actually  biases  the 
search  direction.  This  is  undesirable  since  we  wish  to  use 
the  push-off  factors  as  our  means  of  controlling  the 
search  direction.  A  method  which  avoids  this  bias  in  search 
direction  is  the  circle  as  shown  Figure  2.3  .  The  norm  here 
is 

S-S  <  1  (2.9.1) 

4 .  Simple  Simplex-like  He thod  for  Search  Direc tion 

Vanderplaats  [Bef.  5:  pp. 168-169]  provides  the 
matrix  formulation  which  solves  the  above  sub-optimization 
problem  by  using  the  Toutendijk  method. 


Maximize  P • v  (2.10) 

Subject  to  ; 

A.y  <0  (2.  11) 


20 


y.y  <  1 


(2.  12) 


Wher  e 


(2.  13) 


VG  (X),  0 

2Tg2(X),  q2 


(2.14) 


?G.(X),  q. 

3  °3 

.  VF  (X),  1 

and  where  j  is  the  number  of  active  constraints  (JJAC) 

When  the  solution  to  Eg  (2. 10)  through  (2.12)  is  found,  S  may 
be  normalized  to  some  value  other  than  unity,  but  the  form 
of  the  normalization  is  the  same.  A  solution  to  the  above 
problem  may  be  obtained  by  solving  the  following  system 
derived  from  the  Kuhn-Tucker  conditions  for  that  problem  : 


B  I 


u  >0 
i  - 


v  >0 

i  - 


u  •  v  =  0 


(2.  15) 


(2. 16) 


Wher  e 


(2.  1 7  > 


I  =  Identity  matrix  (2.18) 

C  =  -A-F  (2.  19) 

Above  system  can  be  solved  using  a  complimentary  pivot  algo¬ 
rithm.  Choose  an  initial  basic  solution  to  Eg (2.  15)  is  to 
be 


y  =  c ,  u  =  0  (2.20) 

where  v  is  the  set  of  basic  variables  and  u  is  the  set  of 
nonbasic  variables.  If  all  v^  >  0,  Eg (2.  16)  is  also  satis¬ 
fied  and  problem  is  solved.  If  some  v^  <  0,  the  solution 

procedure  is  as  follows  : 

Let  E-,g  be  the  diagonal  element  of  the  i-th  nonbasic  vari¬ 
able. 

1.  Giver  the  condition  that  some  c  is  less  then  zero, 
we  find  max  (ct/3pi.)  which  is  the  incoming  row  to  the 
basis. 

2.  The  incoming  column  is  changed  to  a  basic 
column,  the  tableau  is  updated  by  a  standard  simplex 
pivot  on  9^  . 

3.  Until  all  c^>  C,  repeat  steps  1.  and  2. 

4.  Tlhen  all  >  0,  the  iteration  is  complete.  The 

value  of  u  is  now  the  desired  solution. 

5.  3y  using  y  =  p-AT'U,  we  get  the  usable-feasible 
search  direction  3  which  is  first  NDV  components  of 

y- 

c •  Initially  lBlS§sible  Designs 

The  method  of  feasible  directions  assumes  that  we 
begin  with  a  feasible  design  and  feasibility  is  maintained 
throughout  the  optimization  process.  If  the  initial  design 
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is  infeasible,  then  a  search  direction  pointing  toward  the 
feasible  region  can  be  found  by  a  simple  modification  to 
direction  finding  problem. 

A  design  sit  nation  can  exist  in  which  the  violated 
constraints  are  strongly  dependent  on  part  of  the  design 
variables,  while  the  objective  function  is  primarily  depen¬ 
dent  on  the  other  design  variables.  This  suggests  a  method 
for  finding  a  search  direction  -which  will  simultaneously 
minimize  the  objective  while  overcoming  the  constraint 
violations.  These  considerations  lead  to  the  following 
statement  of  the  direction  finding  problem  [Ref.  5  : 

pp. 171-172  ]  : 

Maximize  -  VF  (X)  •  S  ♦  (2-21) 

Subject  to  ; 

tfG  (X)  •  S  ♦  0  £  <  0  j  e  J  (2.  221 

j 

S • S  <  1  (2.  23) 


where  J  is  the  set  cf  active  and  violated  constraints,  and 
where  the  scalar  §  in  Eg  (2. 21)  is  a  weighting  factor  deter¬ 
mining  the  relative  importance  of  the  objective  and  the 
constraints.  Usually  a  value  cf  3L  >  10000  will  ensure  that 
the  resulting  S-vector  will  point  toward  the  feasible 
region.  Incorporating  Eg(2.21)  and  Eg(2.22)  into  the  direc¬ 
tion  finding  algorithm  requires  only  that  we  modify  the 
p-vector  given  in  Eg  (2.  24)  and  the  A-matrix  of  Eg  (2. 25). 


P 


-  SF  (X)  ] 


£ 


(2.24) 
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25  (X),  01 

*VX)'  ©2 


(2.25) 


7G.(X),  0  . 

**  3  3 


0  <50 

3  - 


(2.26') 


He  use  the  simple  simplex-like  method  to  fine  the 
search  direction  toward  the  feasible  region. 


C.  CNE-DIHENSIONAL  SEARCH 

1  •  No  Isolated  C enstrai nts 

If  no  constraints  are  violated,  we  find  the  largest 
a*  in  Eg  (1.4)  from  all  possible  values  that  will  minimize 
the  objective  on  S  without  violating  any  constraints,  active 
or  inactive. 

The  procedure  in  HSCCF  is  as  follows  : 

1.  Let  aO,  al,  a2,  a3  be  the  scalar  in  Eg  (1.4)  corre¬ 
sponding  to  pc_nts  X0  ,  XJ ,  X2,  X3,  X4. 

2.  aC  =  0  at  given  point  X  0. 

3.  In  order  to  get  al,  we  can  calculate  the  al  to 
reduce  the  objective  by  at  most  10%  or  to  change  each 
of  the  design  variable  X.  by  at  most  10%. 

4.  Update  the  design  variables  to  XJ  using  E  g  ( 1 . 4)  . 

5.  Evaluate  the  objective  for  XJ ,  ana  check  the  feasi¬ 
bility.  If  one  or  more  constraints  is  violated,  then 
al  is  reduced  to  al/2,  and  we  go  to  step  4. 

6.  In  order  to  estimate  a2,  we  can  use  the  quadratic 
approximation  with  2  points  X,  XJ  and  the  y?. 


a  1  to 


7.  Update  the  design  variables  to  12  by  Eg  ( 1 . 4)  and 
check  the  side  constraints. 

P.  Evaluate  the  objective  and  constraints. 

9.  Now  having  3  a's,  and  values  of  objectives  an1, 
constraints  for  design  variables  10,  XJ ,  X_2  are 

known,  so  by  using  3-point  quadratic  approximation,  a 
value  of  a 3  is  found. 

1C.  Update  the  new  optimal  point  in  search  direction  by 
Eg  (1.4)  . 

11.  Evaluate  the  objective  and  constraints. 

12.  Now  choose  last  3  values,  al,  a2,  a3  and  find  a  new 
a3  using  3-poir.ts  Quadratic  approximation 

13.  Choose  the  a*  among  the  5  points  which  corresponds  to 
the  minimum  objective  function  value  with  no-viclated 
constrain  ts. 

2 •  Cne  or  More  Constraints  Violated 

If  one  or  more  constraints  are  initially  violated,  a 
modified  usable-feasible  direction  is  found.  It  is  then 
necessary  to  find  the  scalar  a*  in  Eg  (1.4)  which  will  mini¬ 
mize  the  maximum  constraint  violation,  using  the  most 
violated  constraint  j,  a  good  initial  estimate  for  a*  is 

“G  .  (X) 

a*  =  - - -  (2.27) 

VG  (X)  •  S 
3 

Since  the  gradients  of  the  violated  constraints  are 
known,  the  scalar  which  is  required  to  obtain  a  feasible 
design  with  respect  to  violated  constraint  in  the  search 
direction,  is  given  to  a  first  approximation  by  Eq(2.27). 

The  more  detail  procedure  in  XSCOP  is  as  follow  ; 

1.  Choose  the  most  violated  constraint  j. 

2.  Calculate  a*  for  violated  constraint  j  usinj 
Eg  (2.27)  . 


3. 


Update  the  design  variables  for  a*  and  check  the  side 
constrain  ts. 

4.  If  one  or  more  violated  constraints  still  exist,  then 
calculate  the  derivative  of  objective,  violated  and 
active  constraints  and  find  a  new  search  direction 
and  then  go  tc  step  1.  Otherwise  proceed  with  the 
optimization  in  the  normal  fashion. 

E.  CONVERGENCE  CRITERIA 

A  desired  property  of  an  algorithm  for  solving  a  nonli¬ 
near  problem  is  that  it  should  generate  a  sequence  of  points 
converging  to  a  global  optimal  point.  In  many  cases, 
however,  we  may  have  to  be  satisfied  with  less  faverable 
outcomes.  In  fact,  as  a  result  of  non-convexity,  problem 
size,  and  other  difficulties,  we  may  stop  the  iterative 
procedure  if  a  point  belongs  to  a  described  set,  which  is 
defined  in  USCOP  as  fellows  ; 

^  Q1  =  U  i  I*0  -  II  <  £*II°I  ) 

2.  Q  =  {XI  |F(r>)  -  F{X$1  <  £.•  I F  (X° )  |  } 

2  * 

In  M SCO F ,  the  algorithm  is  terminated  if  a  point  X_  is 
reached  such  that  X  6  Q,  fl  Qz  .  is  0.03  1  and  is 

approximatly  0.001.  Since  in  engineering  design  problems  it 
is  net  necessary  to  find  solutions  with  more  than  three 
significant  digits. 
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III.  HSCCP  USAGE 


A.  INTRODUCTION 

Since  this  MSCO?  is  written  in  WATERLOO  BASIC  Version 
2.0,  it  is  very  convenient  to  use.  The  user  must  first 
formulate  the  design  problem  with  the  classical  machine 
design  criteria.  Given  the  formulation  of  the  design 
problem  as  a  nonlinear  program,  the  user  then  enters  the 
problem  as  a  part  of  a  BASIC  program.  The  user  defines  the 
objective  function  and  constraint  functions  using  EASIC 
statements.  Other  parameters  are  input  as  data  :  the  number 
of  design  variables  NDV,  the  number  of  inequality 
constraints  NIQC,  variable  bounds  an  initial  design  value 
and  a  print  control  number. 

E.  PRCBIEF1  FOBMOLATICN 

Generally,  the  experienced  design  engineer  will  be  able 
to  choose  the  appropriate  objective  for  optimization 
depending  on  the  requirements  of  the  particular  application. 
The  physical  phenomena  of  significance  should  first  be 
summarized  for  the  device  to  be  designed.  The  appropriate 
objective  can  then  be  selected  and  constraints  can  be 
imposed  cn  the  remaining  phenomena  to  assure  an  acceptable 
design  from  all  standpoints.  However,  the  initial  formula¬ 
tion  for  the  optimization  problem  should  not  be  more  compli¬ 
cated  then  necessary  and  this  often  requires  the  making  of 
some  simplifying  assumptions.  [Ref.  9]. 

After  completing  the  formulation  of  the  design  problem, 
the  design  engineer  should  be  able  to  answer  the  following 
questions  : 

1.  What  are  the  design  variables  ? 


2.  What  is  the  objective  function  ? 

3.  What  are  the  inequality  constraints  ? 

4.  What  are  the  bounds  on  the  variables  ? 

The  engineer  is  then  ready  to  input  the  program  to  the 
MSCOF.  However,  additional  study  and  preparation  of  the 
problem  may  be  useful.  In  particular,  redundant  constraints 
should  be  avoided  if  possible.  MSCOP  will  operate  with 
redundant  constraints  but  it  will  operate  faster  without 
them.  Selection  of  an  initial  design  point  from  which  to 
start  this  program  is  important,  since  it  affects  perform¬ 
ance  and  running  time.  The  user  should  use  any  available 
information  which  gives  a  good  initial  approximation.  If 
side  constraints  exist,  the  user  must  be  sure  the  initial 
values  cf  the  design  variables  do  not  violate  the  side 
constraints.  This  program  will  automatically  handle  an 
initial  design  point  which  is  infeasible  with  respect  to  the 
G(X)  <0  constraints.  However,  if  the  initial  point  does  not 
violate  these  constraints,  convergence  will  likely  be  more 
rapid . 

C.  PECBIEM  SNTEY 

Problem  entry  is  accomplished  by  editing  the  main 
program  directly.  As  an  example,  consider  the  following 
simple  HIP  with  two  design  variables,  and  three  constraint 
funct ions. 

2  2 

Minimize  F(X)  =  X  +  3  X  X  +  2X-X-X+1 

1  12  2  12 

subject  to  ; 

X  +  X  -  3  <  0 
1  2 

1  1 

-  + - 2  <  0 

XX- 
1  2 


28 


X  +  X-  X-  2<0 

1  12- 


X  >0.1 
i  - 

With  the  MECOP  loaded  into  memory  and  listed  on  the  CET, 
modifications  are  made  on  the  program  lines  as  follows  to 
input  this  example  : 

Line  100 

Just  after  the  word  "data”,  three  integers  are  added, 
separated  by  a  comma.  The  first  number  is  11DV  which  is 
the  number  of  design  variables,  the  second  is  MIQC  which 
is  the  number  of  inequality  constraints,  and  the  third  is 
IFET  which  is  print  control  number  (  0  ;  only  final 

results,  1  ;  given  data  and  final  results,  2  ;  given  data 
and  iterative  subcptimal  results) 

for  example  : 

100  data  2,3,2 

Lines  201-220 

Each  line  here  corresponds  to  a  separate  design  variable, 
beginning  with  X  (1)  and  continuing  in  order  to  input 
X  (ND V) .  On  each  line,  three  values  are  separated  by 
commas.  After  the  word  "data",  these  values  are  the 
initial  values  of  the  design  variable,  the  lower  bound  on 
the  variable  and  the  upper  bound  on  the  variable.  Tf  no 
bound  is  to  be  specified,  the  entry  is  filled  by  "no". 

For  the  sample  problem,  the  input  is  : 

2C1  data  3.  ,0.1, no 
202  data  3.  ,  0 . 1,  no 


lines  400  -  450 


These  lines  are  available  for  defining  the  objective 
function.  The  objective  function  East  be  defined  in 
terns  of  subscripted  design  variables  X  ( 1) ,  7(2),  etc. 

Tor  the  sample  problem,  the  input  is  : 

400  f  n_f  =  x  ( 1)  **2  +  x  (1)  *x  (2)  +2.  *x  (2)  **2-x  {  1) -x  (2)  +1 . 

Lines  500-650 

These  lines  are  available  for  defining  the  inequality 
constraint  functions,  which  must  be  expressed  using  the 
format  : 

601  if  i  =  k  ther  fn_g  =  G  (x)  -  b 

i  i 

Tor  the  sample  problem,  the  input  is  : 


OC601  if  i 
00602  if  i 
C0603  if  i 


1  then  fn  g  =  x(1)+x(2)-3. 

2  then  fn_g  =  1 .  /x  (1)  +  1 . /x  (2)  -  2 . 

3  then  fn_g  =  x  { 1)  **2  +  x  (1 )  -x  (2) -2. 


If  there  are  many  constant  values  in  the  constraint  func¬ 
tions,  the  user  may  input  data  for  these  functions  on 
lines  501-600  in  order  to  simplify  their  statements. 
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kt-  IV.  EXAMPLE  PROBLEMS 


A.  DESIGN  0?  CANTILEVERED  BEAM 

1.  Uniform  Cantilevered  Beam 

Assume  a  cantilevered  beam  as  shown  in  Figure  4.1 
must  be  designed.  The  objective  is  to  find  the  minimum 


volume  of  material  which  will  support  the  load  P. 

The  design  variables  are  the  width  B  and  height  H  in 
the  team.  The  design  task  is  as  follows  :  Find  B  and  H  to 

minimize  volume  V  =  B  H  1  (4.  1) 


3  1 


we  wish  to  design  the  team  subject  to  limit  on  bend  in 
stress,  shear  stress,  deflection  and  geometric  conditions 
The  Lending  stress  in  the  beam  must  not  exceed  20,000  psi . 


Me  6  P  1 

<T  =  -  =  -  <  20,000  (4.  2) 

b  I  2  - 

3  H 


The  shear  stress  must  not  exceed  10,000  psi. 


3  P 
2  A 


3  P 
2  B  H 


<  10,000 


(«•  3) 


and  the  deflection  under  the  load  must  not  exceed  1  inch. 


3 

4  F  1 


3 

E  B  H 


<1.0 


(4.  4) 


Additionally,  geometric  limits  are  imposed  on  the  team  size 
0.5  <  B  <  5.0  (4.5) 


1.0  <  H  <  20.0 


(4.6) 


H/b  <  10.  (4.7) 

Now  we  can  input  this  problem  to  MSCOP. 

Input  NDV,  NIQC,  IPKT 

00100  data  2,4,2 

Initial  starting  points 

00210  data  3. 5, 0.5, 5.0 
00220  data  16.0,1.0,20.0 

Evaluation  of  objective 


00400  f r._f  =  tl*x  (1)  *x  (2) 


Evaluation  of  constraints 


00  5  00 

tl 

= 

200 

0  0  c  C  1 

be 

= 

30  . 

e  +6 

0  0  5  0  2 

bp 

= 

10000. 

0  0  5  0  3 

i  r 

i 

=  1 

then 

fn  g  = 

00  5  0  3 

if 

i 

-  0 

then 

fn  g  = 

0050  3 

if 

i 

=  3 

then 

f  n  g  = 

00  5  0  3 

if 

i 

=  4 

then 

t  n  g  = 

6.  *bp*t  1/(20  000  .  *b*h*»  2)  -1  . 
3.  *bp/(  1000  0. *2. *b*h) -1. 

«.  *bjb*tl**3/ (be*b*h**3)  -1. 


TABLE  I 

The  Solution  of  a  Uniform  Cantilevered  Beam 


objective  ;  6664.0 

design  variable  ; 


U 1) 

= 

1.852 

X  (2) 

— 

17.99 

constraint  ; 

gO) 

0.000902 

g(2) 

= 

-0.9549 

g(3) 

= 

-0.0109 

g  («*) 

- 

-0.0286 

As  a  result  of  this  problem  are  in  Table  4.1. 

2 •  VSEiahlg  Cantilevered  Beam 

The  cantilevered  beam  shown  in  Figure  4.2  is  to  be 
designel  for  minimum  material  volume.  The  design  variables 
are  the  width  b  and  height  h  at  each  of  5  segments.  We 
wish  to  design  the  beam  subject  to  limits  on 
stress  (calcula te d  at  left  end  of  each  segment),  deflection 
under  the  load,  and  the  geometric  reguirerae nt  that  the 
height  of  ar.y  segment  does  not  exceed  20  times  the  width. 


33 


where  the  deflection  y  is  defined  as  positive  downward,  y* 
is  the  derivative  of  y  with  respect  to  the  X,  and  1;  is  the 
length  of  of  segment  i.  Young’s  modulus  E  is  the  same  for 
all  segments,  and  the  moment  of  inertia  for  segment  i  is 

3 

t  h 
i  i 

I  =  -  (4.  Ill 

i  12 


The  tending  moment  at  the  left  end  of  segment  i  is  calcu¬ 
lated  as 


1  =  ?  L  +  1-  2I1 

i  i  D=1  i  J 


(4. 12) 


and  the  corresponding  maximum  bending  stress  is 

M  h 


<r 

i 


x  1 

2  I. 


l 

The  design  task  is  now  defined  as 

M 


."liniraize  : 
Subject  to  : 


V  =  ZZ  b 
i=1  i 


h 

i 


1 

i 


O'; 

- - - 1  <  0  i  =  1 ,  .  .  . ,  N 

<7 

ys 

- 1  <  0 

y 

h  -20  b  <  0  i  =  1 ,  .  .  . ,  N 

i  i  - 


(U.  13) 


(4.  14) 


(4.  15) 


(4.  16) 


(4.  17) 


(4.  18) 
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b  >1.0 
i  - 


h  >5.0 

i  - 


i  —  1  /  •  •  •  9 


(4.  19) 


where  ?  is  the  allowable  bending  stress  and  y  is  the  allo¬ 
wable  displacement.  This  is  a  design  problem  in  10  vari¬ 
ables.  There  are  6  nonlinear  constraints  defined  by  Eg  (4.  16) 
and  Eg  (4.  17),  and  5  linear  constraints  defined  bv  Eg  (4. 15), 
and  10  side  constraints  on  the  design  variables  defined  by 
Eg  (4  .  19)  . 

How  we  can  input  this  problem  to  1SC0P. 

Input  NDV,  NIQC,  IPRT 
0010C  data  10,  1  1,2 
Initial  starting  points 


00210 

00220 

00730 

00240 

00250 

00260 

00270 

00280 

00290 

00300 


data  5. 
data  5. 
data  5. 
data  5. 
data  5. 
data  40 
data  40 
data  40 
data  40 
data  40 


1 .  ,no 
1 . ,  no 
1 .  , no 
1 . ,  no 
1 .  ,  no 
, 5. , no 
,5. , no 
, 5. , no 
,5. ,no 
,5., no 


Evaluation  of  objective 

00400  fn  f  =  100.  *  {  x  ( 1 )  * x  (6) 

x(4)  *x(9)  -+  x  (5)  *x(10)  ) 

Evaluation!  of  constraints. 


♦  x  (2)  *x  (7)  +  x  (3)  *x  (3) 


00498  dim  bm7?0)  .bl  (10)  ,sigi  (10)  ,ypb(10)  ,yb(10) 
00500  pet  =  oOOOO. 

00501  be  =  200. e+5 

00502  tl  =  200. 

30503  sigb  =  14000. 

005C4  ytb  =  .5 
00505  si  =  40. 

00506  fer  m  =  1  to  5 

00507  bn  (m)  =  pet*  (t  1  +  s1-id*s1) 

00 508  next  m 

00509  for  e  =  1  to  5 

00510  km  =  m+5  ^ „ 

005  1  1  bi(m)  =  x  (m)  *x  (km)  **  3/1  2  . 

005 1  2  sigi(m)  =  bm  (mj  *x  (km)  /  (2.  *bi  (id)  ) 

00513  next  m 
005  1  4  yzo  =  0 . 

00515  ypzo  =0. 

00516  for  m  =  1  to  5 


005  1  7 
005  1  p 
005  15 
00500 
0052  1 
00522 
0  0550 
00560 
00570 
00560 
00560 
00600 
006  1  0 
00620 
006  30 
00632 
006  34 
00636 


mri  <■>> 

yl  (i)  =  dumm/  (2- *be*bi  (m)  )  +ypzo*sl  +  yzo 

ypzo  =  y pb  (a) 
y  2  o  =  yb(ra) 
next  n 

rein  constraint  function. 


if 

if 

if 

if 

if 

if 

if 

if 

if 

if 

if 


then  tn_g 
then  fn  g 
then  fn~g 
then  fn  g 


then  fn_g 
then  fn  g  = 
then  fn~g  = 

8  then  fn  g  = 

9  then  fn  g  = 

10  then  fn  g  = 

11  then  f n_g  = 


.  1)  /sigb- 1 . 
?(  /sigb-  1 . 
/sig  b-  1 . 
/sigb-  1 . 
,  ,  /sigb-  1 . 
5)  /ybb-1 . 

)  -20  .  * x  1) 

7)  -20  .  *x  2) 
x  (8)  -20.  *x  3) 
x  (9) -2  0  .  *x  (41 
x  (ID)  -20.  *x  (5) 


TABLE  II 

The  Solution  of  a  Variable  Cantilevered  Beam 
objective  ;  62  1  33-  35 


design  variables  constraints 


x  Cl) 

2.994 

G  (1) 

= 

-0. 00219 

X(2) 

= 

2.782 

G(2) 

= 

-0. 00415 

X(J) 

= 

2-526 

G(3) 

= 

-0- 00508 

X<4) 

= 

2-203 

G  (4) 

= 

-0.0040b 

X  (5) 

= 

1.761 

G  ( 5) 

= 

-0-0177 

X(b) 

= 

CD 

CD 

» 

CT 

to 

G(b) 

= 

-0. 4401 

x  (7) 

= 

55.62 

G<7) 

= 

-0-0101 

X  (8) 

= 

50-  56 

G  (8) 

- 

-0.0231 

X<9) 

= 

44-  14 

G  (9) 

= 

0-0000 

X(10) 

= 

35.  19 

G  { 10) 

= 

-0. 0248 

G  (  1 1 ) 

= 

-0. 0278 

ypzo 


I 
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STAPLE  TRUSS 


E. 

Y 


p  =  50000  N 

1 

E  =  200  GPa 

1  ;j^=  +  14000  M/cm 

£.  =  100  cm 

P 


Figure  4.  3  Design  of  a  5-Bar  Truss. 

A  simple  truss  with  5  members  as  shown  in  Figure  4.3  is 
designed  for  the  minimum  volume.  The  design  variables  are 
the  sectional  areas  cf  the  members.  The  constraints  are 
formed  for  the  stresses  of  the  members  not  to  exceed  the 
given  allowable  stress.  The  lower  bound  for  each  design 
variable  is  also  considered.  The  stresses  are  obtained  by 
the  disp la cement  met  hod  of  the  finite  element  analysis.  The 
equation  to  be  solved  is  given  by 

K- u  =  P  ( 4. 20) 


where  K  is  the  stiffness  matrix,  u  is  the  displacement 
vector  and  P  is  the  lead  vector  as  follows  : 


ffi 


A1 


■F* ~ 


2  2 
)  +  (  1  -  v  )  -  1 

2  2  2  4 


2  2 
A1  =  /(  1  +  u  )  +  (  1  +  v  )  -1 

5  W  3  1  2  1  5 


The  design  problem  is  given  by 


minimize  V  = 


Subject  tc 


5 

21  A  1 
i=  1  i  i 


G  = 
i 


I  (Til 


-  1 . 0  <  0  i  =  1 ,  ..  .  ,  5 


(4.26) 


(  *•  27) 


A  >0.1 
i  - 


X  1  ,  m  .  •  ,  5 


(4. 28) 


The  MSCOF  input  for  this  problem  is  given  as  follows  : 

Input  ND7,  NIQC ,  IPRT 

00100  lata  5,5,2 

Initial  starting  point 

00200  data  3.,.i,no 
00202  data  3.,.1,no 
00204  data  3-,.1,no 
00206  data  3.,.1,no 
00208  data  3.,.1,no 

Evaluation  of  objective 

00400  fn  f  =  100  *  (  x  (1)  +  x  (2)  +  x(3)  ♦  sgr(2.)*x(4) 

scr  (2.)  ) 

Evaluation  of  constraints 


0  50  0  dim  vv  (5) 

050  1  te  =  2.e+  7 
0502  tl  =  100. 

0503  sigb  =  14000. 
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0514 
0505 
0506 
0507 
0506 
0509 
0510 
051  1 
0512 
05  13 
05  14 
0515 
0516 
0517 
051  8 
n51  9 
0520 
0  52  1 
0522 
0523 
0590 
0592 
0594 
0  59  5 
0596 
0598 
0600 
0602 
0604 
0606 
0608 
0610 
0650 


cs  =  2.  *sqr  (2  . ) 
ct  =  te/tl 

ui  i  *ct 

k  2  1  =  k  1 2 

k22  =  (x  (2) +x  (5)/cs)  *ct 

k 2 4  =  -x  2) *ct 

k33  =  (x  (3)  +x  (4)/cs)  *ct 

k 34  =  -x  (4  )  *ct/cs 

k  42  =  k2  4 

k43  =  k34 

k  4 4  =  (x  j2)  +x  j  4) /cs)  *ct 

dk 2  =  -  (k  1  2+k22*dk1)  /k24 
dk3  =  -k33*dk2/k34 
pp  =  -50000. 

vv  (1)  =  pp/(k42*dk1+k43*dk3  +  k44*dk2) 
vv  (21  =  a 6  1  ' 1' 

vv  (3)  = 
vv  (4)  =  d 
dll  =  sir 
dl2  =  sqr 
d!3  =  sqr 


<11  =  |^/jk42*dk1+k43*dk3  +  k44*dk2) 

v  (3)  =  dk3*vv  <\) 

(4)  =  dk 2* vv  (1) 

1  =  s  jr  (  (tl  +  v  v  ( 1)  )  *  *  2  +  vv  (2)  **2)  -  tl 

2  =  sqr  (  (t  1  +  v v  (2) -v v  (4)  )  **2  +  { V v  ( 1) 

3  =  sqr j  (tl+ vv (3  )  **2  +  vv  (4)  **2) -tl 

=  r  n  J 


-  vv  (  3)  )  **2) 


rr  —  a  -  j  x  -  r  .  ,  ■  \  - »  /  -  •  ■  v  •  /  -  /  -  * 

hi  =  sqr  (2 . )  *t  1 

d  1 4  =  sqr  {  (hl+vv  (3)  )  **2  +  (hl-vv  (4)  )  **2)  -hi 
dl5  =  sqr  (  (hl  +  vv(1)  j  **2+ (hl  +  vv  (2)  )  **2) -hi 
reir  constraint. 


if  i  = 
if  i  = 
if  i  = 
if  i  = 
if  i  = 
f  nend 


t  hen 
t  hen 
t  hen 
then 
t  hen 


f  n  g 
f  n~g 
f  n  g 
fr.  g 
f  n  g 


te*dl 1/f tl*sigb 
te*dl2/(tl*siab 


te*<313/ 
t  e  *  d  1 4  / 
te*dl5/ 


’  t l*sigb 
|hl*sigb 
hl*sigb 


TABLE  III 

The  Solution  of  a  5-Bar  Truss 


ob  jective 


108.52 


design  variables 

X  (1)  =  0.  1 

X  (2)  =  0.1 
X  (3)  =  3.514 

X  ( 4)  =  4.946 

v  ' 5)  =  n  1 


constraint 

G  ( 1)  =  -  1. 9988 
G  (2)  =  -2.0030 
G  (3)  -  -0.  0030 
G  (4)  =  -0.  120  3 
G  <51  =  -  1.  8797 


V.  SUMMARY  AND  CONCLUSION 


Numerical  optimization  is  a  powerful  technique  for  those 
confronted  with  practical  engineering  design  problems.  It 
is  also  a  useful  tool  for  obtaining  reasonable  solutions  to 
the  classical  engineering  design  problems.  Since  many  engi¬ 
neers  are  now  using  aicrocomputers  for  solving  design  prob¬ 
lems,  the  development  of  microcomputer  software  which  can  be 
easily  used  is  needed. 

In  this  thesis,  an  algorithm  for  constrained  optimiza¬ 
tion  problems  is  programmed  in  standard  BASIC  language 
(NBAS IC  version  2.0)  on  an  IBM  3033.  The  users  can  easily 
convert  this  to  ether  microcomputers. 

MSCOE  (Microcomputer  Software  for  Constrained 
Optimization  Problems)  employs  the  method  of  feasible  direc¬ 
tions  and  specific  modifications  of  a  one-dimensional  search 
for  constrained  optimization.  MSCOP  has  been  validated  by 
tests  on  three  constrained  optimization  problems.  Its 
performance  is  good  and  could  be  made  better  through  refine¬ 
ment  cf  the  algorithm. 

Since  microcomputers  are  available  with  reasonable 
memory  size  and  computational  speed,  their  capabilities  will 
continue  to  improve  as  more  engineering  software  becomes 
available.  MSCOP  is  considered  to  be  a  first  step  toward 
more  widespread  use  cf  optimization  techniques  on  microcom¬ 
puters. 
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APPENDIX  A 

MSCOP  PBOGEAM  LISTING 


00  10 
0020 
0021 
0020 
004  0 
0050 
0060 
0070 
0030 
0090 
0100 
0110 
01  1  5 
0120 
0125 
0130 
013^ 
0140 


option  base  1 

4  4  _  _  A  I 


d  lit  x(2lj  .xO  (2  1)  ,  gcv  (5 1)  #ngcv  (51)  ,  df  (21)  ,  dg  (5 1 , 2  1) 
thta  (5l)  ,  wrky  (51. 51) . _  .  .  _ 


dim  a(51'21|fb(5j.51  ,  p  (21),  y (21)  ,  s  12  11  ,  u  ( 51) c  (5  1) 
dim  lwrk  (5  1 )  ,  g  vrk  (5 1 )  ,  wrk  1  (51)  ,  wrx2  ( o  1}  ,  wrk 3  { 5  1) 
dim  wrku  (5  1)  -  wrkl  (5  1  ,  lo wb  (2 1 )  ,  upr b  (2  1)  ,  lo %  (6)  ,  up  *  i 
rei  input  data 


(6) 


gosuk  10000 

rem  input  number  of  design  variables  and  constraints. 


read  ndv ,  n  igc ,  iprt 


0150 


bnlo  else  lowb(i)  = 
bnup  else  uprb  (i)  = 


0160 
0200 
0  210 
0360 
0370 
0375 
0380 
0390 
0400 
04  10 
0420 
0430 
0440 
0450 
0480 
0490 
0500 
0510 
0520 
053  0 


data  2,4,2 
for  i  =  1  to  ndv 

rem  input  initial  value  of  design  variables 
read  x  (i) 
xO  (i)  =  x  (i) 
if  mqc  =  0  then  160 
read  Io$,up$ 

if  loS  =  ’no'  then  lowb  (i) 
value  (lo$) 

if  upy  =  ’no’  then  uprb  (i) 
value  (up$) 
next  i 

data  3.5,0.5,10. 
data  16. ,1.0,20. 

rem  evalute  the  objective-function 
obj  =  fn  f  (x) 
itri  =  1“ 

rem  objective  function, 
def  fn  f  (x) 

f n_f  =  200. *x  (1) *x  (2) 
fnend 

rem  evaluate  the  constraints 
for  i  =  1  to  niac 

gcv  (i)  =  f  n_g  (x , i) 
next  l 

rem  constraint  functions 
def  fn  q  (x,i) 


tl  =  2TJ0. 
be  =  30.e+6 


bp  =  10000. 

if  i  =  1  then  fn_g  = 


0540 

0550 


if  i  =  2  then  fn_g  = 
if  i  =  3  then  fn_g  = 


0560 
0650 
0700 
07  10 
0720 
0730 
0740 
0750 
07  60 
0770 
0780 
0  78  5 
0790 
0800 


6.  *bp*tl)  / 

) .* x  ( 1) 


*x  (1)  )-1. 


if  i  =  4  then  fn_g  =  x(2)/ 
fnend 

rem  initial  counting  number  input 
ical  =  1 

if  ical  >  3  then  stop 

rem  call  the  optimization  code. 

gosub  2000 

rem  print  results. 

rem 

rem  re-counting  number  input, 
ical  =  ical+1 
if  ical  =  3  then  850 
rem  10$  reduce  the  design  variables, 
for  i  =  1  to  ndv 


20  000.*x(1)*x(2)**2)  -1  . 

3.  *bp)  /  (29000.  *x  (1)  *x  (2)  )  -  1 . 

4.  *bp*tl**3)  / 

(be*x  J1)  *x  (2)**3)-1. 


C  3  1  0 
0320 
08  30 
0  3  4  0 
0850 
0860 
0370 
0880 
0890 
0900 
2000 
200  1 
2002 
200  3 
200*+ 
2003 
2010 
2020 

2030 

2040 

2050 
2060 
2070 
2080 
2090 
2100 
2110 
2120 
2130 
2140 
2150 
2160 
2170 
2180 
2190 
2200 
2210 
2220 
2230 
2240 
2250 
2260 
2270 
2230 
2290 
2300 
2310 
2320 
2330 
234  0 
2350 
2360 
2370 
2380 
2390 
2400 
2410 
2420 
2430 
2440 
2450 
2460 
2470 
2480 
2490 
2500 
3000 


x  (i)  =  0.9*x  (l) 
xO(i)  =  x  (i) 
next  l 
goto  720 

rem  10"  increase  design  variables, 
for  i  =  1  to  ndv 
x  (i)  =  1 .  1*x  (i) 
xd.(i)  =  x{i) 


viola  te 


for  i  =  1  to  ndv 
x  (i)  =  1 .  1*x  (i) 
x0.(i)  =  x{i) 
next  l 
goto  720 

rem  calculate  the  obj.  constraint  fen. 
obj  =  fn  f  (x) 
for  i  =  T  to  nige 
gcv(i)  =  f  n_g  (x,  i) 
next  l 
itrg  =  1 
itrg  =  itrg+1 

rem  calculate  the  number  of  active  and  violate 
constraints, 
gosub  3500 

rem  calculate  the  gradient  of  objective  and 
active  or  violated  constraints, 
gosub  3800 

if  nave  =  0  then  2190 
gosub  3900 

rem  calculate  the  push-off  factors 
gosub  4000 

rem  making  the  matrix  c 
rem  normalized  the  df  (i) 
gosub  4100 

rem  normalized  the  DG{i) 
gosub  4200 

if  nvc  >  0  then  gosub  4400  else  gosub  4600 
rem  calaulate  the  usable-feasible  direction  s  (i) 
gosub  5000 
goto  2230 

rem  normalize  the  df(i) 
for  i  =  1  to  ndv 
s  (i)  =  -(df  (i)  ) 
next  i 

rem  normalize  the  s  (i) 
gosub  5700 

rem  one-dimensicnal  search 

if  nvc  =  0  then  gosub  6000  else  gosub  9000 

rem  update  x  for  alph 

gosub  7000 

gosub  7100 

rem  calculate  new  point  value, 
nobj  =  fn_f (x) 
rem  convergence  test 
gosub  6780 

if  walp  <=  accx  and  delf  <=  dabf  then  2470 
itn  =  itri  +  1 

if  itri  >  mxit  then  print  'check  the  problem' 


df  (i) 


else  gosub  9000 


if  itri  >  mxit  then  print  'check  the 
obj  =  nobj 
for  i  =  1  to  ndv 
x0.(i)  =  x(i) 
next  i 

for  i  =  1  to  nice 
gev  (i)  =  f  n_g  (x,  i) 
next  i 

if  iprt  =  2  then  2460 
gosub  9200 
goto  2010 

rem  print  final  results 
print  ******  final  results  *****  ' 
gosub  9200 
return 

rem  initialize  the  integer  working  array 


to  nice 
=  f  n_g  (x,  i) 

2  then  2460 


30  0  5 
3010 
3015 
3020 
3050 
3055 
3060 
3065 
3070 
3  10C 
3105 
3110 
3115 
3120 
3150 
3155 
3160 
3165 
3170 
3200 
320  5 
3210 
3215 
3220 
3250 
3255 
3260 
3265 
3270 
3275 
3280 
3300 
3305 
3310 
3315 
332  0 
3350 
3353 
3356 
3359 
3362 
3365 
3368 
3  37  1 
3374 
3377 
3330 
3383 
3400 
34  Oc 
3410 

34  1  c 
3420 
3425 
3430 
3450 
3455 
3460 
3465 
3470 
3475 
3480 
3500 

350  2 
3  50  4 
3510 
3520 

35  3  0 


for  i  =  1  to  niqm 
iwrk  (i)  =0 
next  i 
return 

rem  initialize  the  integer  working  array 
for  i  =  1  to  niqm 
ivrk(i)  =  0 
next  i 
return 

rem  initialize  the  one-dimension  working  array 
for  i  =  1  to  niqm 
wrkl  (i)  =  0  . 
next  i 
return 

rem  initialize  the  one-di mer.sion  working  array 
for  i  =  1  to  niqm 
wr k2  (i)  =  0  . 

next  i 
return 

rem  initialize  the  one-dimension  working  array 
for  i  =  1  to  niac 
vrk3  (i)  =  gcv  (i) 
next  i 
return 

rem  initialize  the  two-dimension  working  array 
for  i  =  1  to  niqm 
for  j  =  1  t c  ndvm 
wrky(i,j)  =  0. 
next  g 
next  i 
return 

rem  initialize  the  derivative  of  objective  DF  (i) 
for  i  =  1  to  ndvm 
df  (i)  =0. 
next  i 
return 

rem  initialize  the  a  (i  ,  j) ,  p  (i)  ,y  (i)  ,c  (ii 
for  i  =  1  to  ndvm 


i&’lVl  tc 
au,i)  =  0. 


a  (1 

next  g 
next  i 


for  j  =  1  to  niqm 
c(  j)  =  0. 


next  g 
return 

rem  initialize  the  derivative  of  constraints  DG(ifj) 
for  i  =  1  to  niqm 
for  j  =  1  tc  ndvm 
dg  (i,  j)  =  0. 
next  g 
next  i 
return 

rem  initialize  the  bfi,j) 
for  i  =  1  to  niqm 
for  j  =  1  to  niqm 
b  (i #  j)  =0. 
next  g 
next  i 
return 

rem  Calculate  the  number  of  active  and  violate 
constraints, 
gosub  3000 
gosub  3100 
nac  =  0 
n  vc  =  0 

for  i  =  1  to  niqc 


3540 
3550 
3560 
3570 
3530 
3530 
36  10 
3620 
3630 
'*640 
3650 
3660 
3670 
3680 
3690 
3700 
3710 
3720 
3730 
3740 
3750 
3790 
3300 
3805 

38  10 
3815 
3820 
3825 
3830 
3335 
3840 
3850 
3860 
390  0 
390  1 

39  10 
3915 
3920 
3925 
3930 
3935 
3940 
3945 
3950 
3955 
3960 

3  966 
4000 

40  1  0 
40  2  0 
4030 
4040 

4  09  0 
4100 
4  102 
4105 
4110 
411  5 
4120 
4125 
4127 
4130 
4135 
4  140 
4145 
4200 
4202 
4205 
4210 


if  gcv  (i)  >=  vcc  then  3580 
if  gcv  (i)  <  acc  then  3590 

nac  =  nac  +  1 
goto  3590 
nvc  =  nvc+1 
next  i 

nave  =  nac+nvc 
if  nave  =  0  then  3790 
ii  =  1 

for  x  =  1  tc  mqc 

if  gcv(i)  >=  vcc  then  3720 
if  gcv  (i)  <  acc  then  3750 
lvrk  f  nvc  +  iil  =  i 
wrkl  (nvc+ii)  =  gcv(i) 
ii  =  ii+1 
goto  3750 
lvrk  ( j  j)  =  i 
wrk  1  (33'  =  9CV  (-1-) 


wrk  1  (g  j)  =  gcv  (i) 

33.=  33+1 
next  i 
return 

rem  calculate  the  gradient  of  f  (x) 

qosub  3300 

for  i  =  1  to  ndv 

dxi  =  fdm*abs  (x  (i) ) 

if  dxi  <=  mfds  then  dxi  =  mfds 
x  ( i )  =  x  (i)  +dxi 
dobj  =  fn  f  |x) 
df  (i)  =  (doti-ob  j) /dxi 
x  (l)  =  xO  (i) 
next  l 
return 

rem  calculate  the  DG(i,j) 
qosub  3400 
for  i  =  1  to  ndv 
dxi  =  fdm*xfi) 

if  dxi  <  mfds  then  dxi  =  mfds 
x  (i)  =  x  fi)  +dxi 
for  j  =  1  tc  nave 
k  =  iwrk  (j) 
dcon  =  fn  g(x,k) 
ig  ( j,i)  =  (acon-wrkl  ( j)  )  /dxi 
next  a 

x  (i)  =  xO  (i) 
next  i 
return 

rei,  calcilate  the  push-off  factor 
for  i  =  1  to  nave 

thta(i)  =  thtO*  (1. -wrkl  (i)  /acc)  **2 
if  tnta  (i)  >  thta  then  thta  (i)  =  t 
next  i 
return 

rem  normalize  the  DF(i) 
gosub  3200 
fsg  =  0. 

for  i  =  1  to  ndv 

fsg  =  fsq  +  df(i)**2 
next  x 

fsg  =  sgr(fsc) 

if  fsc  =  0.  then  fsg  =  zro 

for  i  =  1  to  ndv 

vrk3(i)  =  ( l./f  sg)  *df  (i) 
next  i 
return 

rem  normalize  the  DG(i) 

gosub  3250 

for  i  =  1  to  nave 


qsq  = 


4215 
4220 
42  2  c 
4230 
4232 
4235 
4240 
4245 
4250 
4255 
4400 
4405 
4410 
4420 
4430 
4440 
4450 
4450 
4470 
4480 
44  9  0 
4500 
4510 
4520 
4530 
4540 
4550 
4560 
4570 
4580 
4596 
4590 
4600 
4605 
46  10 
4620 
4630 
4640 
4650 
4660 
4670 
4680 
4690 
4700 
4710 
4720 
4730 
474  0 
4750 
4760 
4770 
5000 
5002 
5005 
5010 
5040 
5050 
5060 
6  07  0 
5080 
5090 
5100 
5110 
5120 
5130 
5140 
5150 
5160 
5170 


for  j  =  1  tc  ndv, 

gsq  =  gsq+dg  (i,  j)  **2 
next  o 

if<igsqS=r^?’"^en  gsq  =  zro 
for  j  =  1  to  ndv 

wrky  (i,  j )  =  ( 1. /gsq)  *dj  (i,  g) 
next  j 
next  i 

return  .  ,  , 

rem  exist  the  violate  constraints 

qosub  3350 

for  i  =  1  to  nave 

for  i  =  1  tc  ndv  ,  . 

a  (i  /  3)  =  wrky(i,g) 
next  j 

a  (i ,na v  +  1)  =  thta(i) 

next  i 

for  i  =  1  to  ndv 
p  (il  =  -wrk3(i) 
next  i 

p  (ndv+1)  =  phid 
for  i  =  1  to  nave 

RrY-  1  tc  ndv  +  1 
xx  =  a  (i  ,  j)  *p  (j) 

F=  yy+xx 

j 

c  (i)  =  (-i.)*yy 

next  i 
ndt  =  nave 

return  .  t .  .  .  . 

rem  only  exist  active  constraints 
osub  3350 
or  i  *  1  to  nave 


:or 


i  =  1  tc  ndv  . 
ati,j)  =  wrky  (i,g) 

^i,  n<?v  + 1 )  =  thta(i) 


next 
a 

next  i 
for  i  =  1  to  ndv 

a  (navc  +  1,  j)  =  wrk3  (g) 
next  3  , 

a (navc+1 ,ndv+ 1)  =  1. 
p  (ndv  +  1)  =  1. 
for  i  =  1  to  navc  +  1 

cc  =  a  (i,  ndv+1)  *p  (ndv  +  1) 
c(i)  =  (-1.)*cc 
next  i 

ndh  =  navc+1 
return 

rem  calculate  the  usable-feasible  direction 

gosub  3000 

gosub  3250 

gosub  3450 

for  i  =  1  to  ndt 

for  j  =  1  to  ndv+1 
wrky  (j,i)  =  a  (i,  j) 
next  j 
next  i 

for  i  =  1  to  ndt 
for  i  -  1  to  ndb 

ff  =  0.  a  , 

for  k  =  1  to  ndv+1 

tf  =  a  (i.  k)  *vrky  (k,  g) 
ff  =  fr+tr 
next  *  „  ^  , 

b  (i ,  j)  =  (-1-)*ff 
next  g 
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5130 
5190 
520  0 
5210 
5220 

52  3  0 
5240 
5250 
5260 
5270 
5280 
5290 
5300 

53  10 
5320 
5330 
5340 
5350 
5360 
537  0 
5380 
53  35 
5400 
5405 
5410 
5420 
5430 
544  0 
5460 
5470 
5480 
5490 
5500 
5505 
5510 
5520 
5525 
5530 
554  0 
5550 
5560 
5565 
5570 
5580 
5600 
5610 
c>62  0 
5630 
5640 
5650 
5660 
56^0 
5680 
5690 
5700 
57  10 
5720 
5730 
5740 
57  50 
5755 
5760 
5770 
5780 
5820 
6000 
6005 
601  0 
6  015 


next  l 

iter  =  0 

nmax  =  5*ndb 

re®  begin  iteration 

iter  =  iter+1 

cbmx  =  0 . 

ichk  =  0 

for  i  =  1  to  ndt 
ci  =  c  { i) 
bii  =  b  (l.i) 
if  bii  =  0.  then  5340 
if  ci  >  0.  then  5340 
cb  =  ci/tii 

if  cb  <=  cbmx  then  5340 
ichk  =  i 
ctmx  =  cb 
next  i 

if  cbmx  <  zro  or  iter  >  nmax  then  5550 
if  ichk  =  0  then  5550 
jj  =  iwrk(ichk) 

if  =  0  then  iwrk  (ichk)  =  ichk  else  iv rk  (i  chk)  =  0 
if  b  (ichk,  ichk)  =  0.  then  b  (ichk,  ichk)  =  zro 
tb  =  l./fc  (ichk, ichk) 
if  bb  =  0.  then  bb  =  zro 
for  i  =  1  to  ndb 

b  (ichk,i)  =  bb*b(ichk,i) 
next  i 

c(ichk)  =  cbmx 
for  i  =  1  to  ndb 

if  i  =  ichk  then  5530 
bbi  =  t(i,ichk) 
b  (i,  ichk)  =  0  . 
for  j  =  1  to  ndb 
if  i  =  ichk  then  5520 

b  (i ,  j)  =  b  (i  ,  j)  -bbi*b  (ichk,  j) 
next  j 

c(i)  =  c  (i)  -bbi*cbnix 
next  i 
goto  5220 
ner  =  0 

for  i  =  1  to  ndt 
u (i)  =  0. 
j  =  iwrk(i) 

if  i  >  0  then  u  (i)  =  c  (j) 
next  i 

for  i  =  1  to  ndt 
ff  =  0. 

for  i  =  1  tc  ndb 

ff  .=  ff+vrky  (i,  j)  *u(  j) 
next  i 

|it)  :  ?fr££ 

next  l 
return 

rem  normalized  the  s(i) 
ssq  =  0. 

for  i  =  1  to  ndv 

ssq  =  ssq  +  s  (i)  **  2 
next  l 

ssq  =  sqr  (ssq) 

if  fslp  =  0.  then  fslp  =  zro 
for  i  =  1  to  ndv 

s  (i)  =  ( 1 .  /ssq)  *s  (i) 
next  i 
return 

rem  or.e-dimensicnal  search  for  initial  feasible  ocir.t 

rem  calculate  fcr  slope  of  f(x) 

fslp  =0. 

for  i  =  1  to  ndv 


6020 
6  0  2  5 
6005 
604  0 
6045 
6  0  5  C 
6055 
606  0 
6  06  5 
6  06  7 
6070 
6  07  6 
6  076 
6  0  0  0 
60  3  5 
6  0  Q  0 
6  0  9  c 
6100 
6  10  5 
6110 
6115 
6120 
6125 
6130 
6135 
6140 
6145 
6150 
6155 
6160 
6165 
6170 
6  1  7  c 
6130 
6185 
6190 
6  1 9  c 

6200 
6205 
6210 
6215 
6220 
6225 
6230 
6215 
6237 
6240 
6245 
6250 
6252 
6255 
6260 
6265 
6270 
6275 
6280 
6  28  5 
6290 
62  9^ 
6300 
6305 
6310 
6315 


6  32  1 
6325 
6  3  26 


fslp  =  fslp+df  (i)  *s  (i) 
next  i 

rea  idenfy  the  initial  point, 
a  lew  =  0. 
flew  =  obj 
for  i  =  1" to  nice 
wrkl(i)  =  gcv(i) 
next  i 

rent  find  a  1st  ;  the  1st  mid-point, 
if  fslp  =  0.  then  fslp  =  zro 
alst  =  abo  j*f  low/abs  (  fslp) 
for  i  =  1  to  ndv 

if  s(i)  =  0.  then  s  (i)  =  zro 
wain  =  alpx*x  (i)  /abs  (s  (i)  ) 
if  walp  >  alst  then  6  095 
alst  =  walp 

next  i 

rem  update  x  for  alst. 
alph  =  alst 
gosub  ^000 
gosub  7100 

rem  calculate  the  fist  and  wrkl(i) 
f  1st  =  f  n_f  (x) 
for  i  =  1  to  nice 

wrkl  (i)  =  f n_g  (x,i) 
next  i 

rem  check  the  feasibility, 
ncvl  =  0 

for  i  =  1  to  niqc 

if  wrkl  (i)  <  vcc  then  6  170 
ncvl  =  ncvl+1 

next  i 

if  ncvl  =  0  then  6200 
alst  =  0 . 5*a1st 
goto  6105 

rem  find  a2nd  ;  the  2nd  mid-point, 
rem  2-points  quadratic  fit  interpolation 
for  minimum  f  (alpa) . 
aO  =  flow 
al  =  fslp 

a2  =  (f  1st-a1*a  Ist-aO) /(a1st**2) 
if  a2  <=  0.  then  a2  =  zro 
a2nd  =  -al/  (2.*a2) 

rem  2-points  linear  interpolation  for  g  (alpa)  =9, 
for  i  =  1  to  niqc 
a9  =  wrkl  (i) 

if  alst  =  0.  then  alst  =  zro 


al  =  (wrkl  (i) -aO) /a  1st 
if  al  <=  0.  then  al  =  5 


walp  =  -aO/al 

if  walp  <=  0.  then  walp  =  1000. 
if  walp  >=  a2nd  then  6265 
a2na  =  walp 

next  i 

rem  update  x  for  a2nd. 
alph  =  a2nd 
gosub  7000 
gosub  7100 

rem  calculate  f2nd  and  wrk2  (i) 
f  2nd  =  fn  f  (x) 
for  i  =  1~to  nice 

wrk2  (i)  =  f n_g  (xri) 
next  i 

rem  find  final  point  aupr  by  using 
3-points  quadratic  fit. 
fl  =  flow 
alpl  =  alow 
f2  =  fist 
alp2  =  alst 
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6330 

63  3  1 
6335 
6340 

6342 

6343 
634  7 
6350 
6355 
636  0 
6365 
63^0 
6  375 

6376 

6377 
6  3  8  C 
6  38 c 
6390 
63S5 
6400 
6405 
6410 

64  1  5 
6420 
6425 
6430 
6435 
6440 
6445 
6450 
6455 
6460 
6465 
6470 
6475 
6480 
6485 
6490 
6495 
6500 
6505 
6510 

65  15 
6520 
6525 
6530 
6535 
6540 
6  54  5 
6550 
6555 
6  56  0 
6565 
6567 
6  56  9 
6  57  1 
6573 
6575 
6577 
6575 
6600 
669  3 

660  5 

66  10 
66  1c 
662  0 
66  30 


£3  =  f 2nd 
alp3  =  a2nd 

rem  3-ooints  quadratic  fit  interpolation 
gosub  5600 

if  a2  =  0.  then  a2  =  zro 
a3rd  =  -a  1  /  f  2.  *a2) 
if  a3rd  <=  0.  then  a3rd  =  1900. 
for  i  =  1  to  niqc 
fl  =  wrklfi* 
f 2  =  wrkl  (i 
f  3  =  wrk2]i 
gosub  660C 
gosub  6630 

if  alps  >  a3rd  then  6380 
a3rd  =  alps 

next  i 

rem  update  x  fcr  aupr 
alph  =  a  3rd 
gosub  "’OOO 
gosub  7100 

rem  calculate  the  fupr  and  wrku(i) 
fupr  =  fn_f(x) 
for  i  =  1  to  nigc 

wrku  (i)  =  f  n  g  (x,i) 
next  i 

rem  find  4th  nev  point, 
fl  =  fist 
f  2  =  f 2nd 
f 3  =  f 3r d 
alpl  =  a  1st 
alp2  =  a2nd 
alp3  =  a3rd 

rem  3-points  quadratic  fit. 
gosub  5600 

if  a2  =  0.  then  a2  =  zro 
aupr  =  -a1/(2.*a2) 
for  i  =  1  to  nigc 
f 1  =  wrkl  (i) 
f  2  =  wr  k2  (i 
f3  =  vrk3  (i 
alpl  =  alst 
alp2  =  a2nd 
alp3  =  a3rd 
gosub  6600 
gosub  6630 

if  alps  >  aupr  then  6540 
aupr  =  alps 
next  i 

rem  update  x  for  aupr 
alph  =  aupr 
gosub  7000 
gosub  7100 

rem  evaluate  fupr  and  wrku(i1 
fupr  =  f  n_f  (x) 
for  i  =  1  to  nigc 

wrku  (i)  =  f n_g  (x, i) 
next  i 

rem  find  optimum  alpa  for  not  violating 

gosub  14300 

return 

re®  3-points  quadratic  fit. 

if  alDi  =  alp2  cr  alp2  =  alD3  or  alpl  = 

then  return 

f  3-f  1)  /  (alp3-alp11  - 
2-f  1)  /  (alp 2- alp  1}  )/Jalp3-alp2) 
(f2-f1)/(alp2-alp1)-a2*(alp1  +  alp2) 
1-a 1 *a lp  1-a2*alp  1**2 

return 

rem  zero  of  polynomial  fcr  g(alpa) 


cons  train  ts . 

alp  3 
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66  35 

664  0 
6642 
6645 

665  0 
6655 
6660 
6665 
6670 
6675 
6680 
6  68  5 
6690 
6695 
6700 
6805 
6710 
6712 

67  15 
6720 
6780 
6790 
6795 
6800 
6915 
6816 
6820 
6330 
6850 
6855 
6860 
6870 
6380 
6890 
6910 
6990 
70  0  0 
7010 
7020 
7030 
7040 
7100 
711  0 
7120 
7130 
7140 
7160 
80  0  0 
8010 
8020 
8030 
8040 
8050 
8060 
8070 
8080 
8090 
8100 
8110 
8120 
9130 
8140 
8150 
8160 
8170 
8180 
8190 
8200 
8210 


di  =  a 1**2-4. *dz*au 
if  dd  <  0-  then  6715 
if  a2  <=  0.  then  a2  =  zro 
if  a2  =  0.  then  a2  =  zro 

alpb  =  {-al  +sgr  (dd)  )  /  (2.  *a2) 
alpc  =  (-al-sgr  jdd{  )/ (2.  *a2) 
if  alpb  <=  0  and  alpc  <=  0.  then  6715 
if  alpb  >=  0.  and  alpc  >=  0.  then  6695 
if  alpb  >=  0.  and  alpc  <  0.  then  6685 
alps  =  alpc 


goto  6720 

alps  =  alpb 
goto  6720 

if  alpb  >=  alpc  then  6710 
alps  =  alpb 
goto  6720 

alps  =  alpc 
goto  6720 
alps  =  1000. 
return 

rem  update  aboj  and  alpx 

delf  =  abs  (obi-nob j) 

diff  =  abs  (delf/objl 

abcj  =  (abo j  +  dif f) /2. 

valp  =  0. 

welx  =  0. 

for  i  =  1  to  ndv 

delx  =  abs  (xO  (i) -x  (i)  ) 
difx  =  abs (delx/xO  (i) } 
if  delx  >=  welx  then  welx  =  delx 
if  difx  <=  walp  then  6880 
walp  =  difx 

next  i 

alpx  =  (alpx+ walp) /2. 
datf  =  accf*abs  (ofc j) 
return 

rem  update  the  x  (l) 
for  i  -  1  to  ndv 

x (i)  =  xO (i) +alph*s  (i) 
next  i 
return 

rem  check  the  side-constraints, 
for  i  =  1  to  ndv 

if  x(i)  <  lcwb(i)  then  x  (i)  =  lowb  (i) 
if  x(i)  >  uprb  (i)  then  x  (i)  =  uprb  (i) 
next  i 
return 

rem  estimate  the  alpa 
fstr  =  flow 
alpa  =  alow 
nvcl  =  0 

for  i  =  1  to  nigc 

if  wrkl  (i)  <  vcc  then  8070 
nvcl  =  nvcl+1 


next  l 
if  nvcl  > 
if  fist  > 


0  then  8120 
fstr  then  8120 


alpa  =  alst 
fstr  =  fist 
nvcl  =  0 

for  i  =  1  to  nigc 

if  wrk2(i)  <  vcc  then  8160 
nvcl  =  nvcl+1 
next  i  „ 

if  nvcl  >  0  then  8210 
if  f 2nd  >  fstr  then  8210 
alpa  =  a2ni 
fstr  =  f2nd 
nvcl  =  0 


8220 
°  2  3  0 
8240 

8  2  5G 
8260 
827  0 
82  8  0 
82^0 
8300 
83  10 
«32  0 
8330 
8340 
8350 
8360 
33  7  0 
83  3  0 
9390 
8400 
9000 

9002 

9004 

9006 

9098 

9010 

9012 
9014 
9016 

9013 

9  0  20 
9022 
7024 
902  6 
902  7 
9028 
90  30 
7032 
90  34 
9036 
9038 
9040 
7042 
9044 
9046 
9048 
9050 
9052 
9054 
9056 
9058 
9060 
9062 
9064 
9066 
9068 
9070 
9072 
9074 
9076 
9078 
9030 
9200 
9205 
9210 
9215 
92  2  0 
9225 
9230 


for  i  =  1  to  rice 

if  wrk3  (i)  <  7cc  then  3250 

nvcl  =  nvcl+1 
next  i 

i c  nvcl  >  0  ther  8300 
if  f3rd  >  fstr  then  9300 

alra  =  a3rd 

fstr  =  f3rd 

nvcl  =  0 

for  i  =  1  to  nice 

if  vrku{i)  <  vcc  then  8340 
nvcl  =  nvc 1 + 1 
next  i 

if  nvcl  >  0  then  8390 
if  f upr  >  fstr  then  8390 

alia  =  aupr 

fstr  =  fupr 

alph  =  alpa 
return 

rent  one-dimensicnal  search  for  initial 
infeasible  point, 
l  i  =  1 

gevm  =  wrkl  ( 1) 
tor  i  =  1  to  nave 
if  wrkl  (i)  <=  gevm  then  9014 
ii  =  l 

qcvm  =  wrkl  (i) 
next  i 

rein  calculate  the  slope  of  badly  violation. 

gslp . =  0 . 

for  i  =  1  to  ndv 


£s}p  =  gslp  +  dg  (ii, i)  *s(i) 


rem  calculate  the  alph. 
n  gslp  =  0.  then  gslp  =  zro 
alph  =  -gcvm/gslp 
rem  update  X  for  alph. 
gosub  7000 
gosut  7100 

rem  evalute  the  objective  and  constraint. 

obg  =  fn  f  (x) 

for  i  =  7  to  nige 

gcv(i)  =  f n_g  (x,  i) 
next  i 

rem  calculate  the  NVC. 
gosub  3500 

if  nvc  =  0  then  return 
rem  update  initial  value, 
for  i  =  1  to  ndv 
x0.(i)  =  x  (i) 
next  i 

go«=uba3800ate  ^  ^  /(^9  (i/ j)  and  push-off  factor. 

gosub  3900 
gosub  4000 

rem  normalize  the  d f  fi)  ,  dg  [i,  i) 
gosub  4100  v  ' '  y  1  ' 

gosub  4200 

rem  find  the  search  direction, 
gosub  4400 
gosub  5000 
goto  9000 

rem  print  the  results 
print  ,f 

print  1  ***********  data  ***********1 
print  * ' 

print  ’The  number  of  design  variables  =  *  r.i 

print  jThe  number  of  ineguality  constraints  = 


q235  print  ’The  objective  value  =  *,obj 
9240  print*' 

924  5  print  ******  design  variables  ****** 

9250  tor  i  =  1  to  ndv 

9255  print  ’x(’;i;')  =  '  ,x  (i) 

9260  next  i 
9265  print  ’ ' 

9270  print  ’the  number  of  active  constraints  =  *;nac 
92 7 5  print  ” 

9280  print  ’the  number  of  vionate  constraints  =  ’  ;nvc 
9285  print  ’’ 

9290  print  *****  constraint  value  ***** 

9295  print  •’ 

9900  for  i  =  1  to  nice 

9305  print  •  g  ( *  ;  i ;  • )  =  *  ;gcv  (i) 

9310  next  i 
9315  return 

9500  rem  default  number 

9510  mxit  =  50  !  maximum  iteration  number 

952C  fdm  =  .01  !  finite  difference  step 

9530  mfds  =  .001  !  maximum  absolute  finite  difference  stec 
9540  vcc  =  .004  !  violated  constraint  criteria  (thickness) 

°550  acc  =  -.  1  !  active  constraints  criteria  (thickness) 

9560  thtO  =  1.  !  push-off  factor  multiplier  (theta  zero) 

9570  thtm  =  50.  !  maximum  value  of  push-off  factor 

9580  phid  =  100000.  !  weighting-factor  used  in  direction 

when  infeasible 

9590  accf  =  .001  !  absolute  convergence  criteria 

9600  accx  =  0.001  !  absolute  convergence  criteria.* 

9610  zro  =  .0001  !  defined  zero 

9620  espl  =  .005  !  used  to  prevent  division  by  zero 

9630  bnlo  =  -l.e+70  !  the  value  of  low  boundary 

964  0  bnup  =  l.e+70  *.  the  value  of  upper  boundary 

9650  dalp  =  .01  !  steo  size  of  alpa  in  one-dimensional 

search 

9660  abej  =0.1  !  step  size  for  reduce  objective 

9670  alpx  =  .1  !  reduce  the  design  variable  factor 

9680  ndvm  =  21  !  the  number  cf  maximum  design  variable 
9690  nigm  =51  !  the  number  of  maximum  inegualitv 

constraints 

9700  return 
9800  end 
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